{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "debe61b0",
   "metadata": {},
   "source": [
    "# Spatial Solver Tutorial"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "37306187",
   "metadata": {},
   "source": [
    "`spatialsolver` is a pyne module that contains seven neutron transport equation solvers.\n",
    "The neutron transport equation is a balance statement that conserves neutrons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "f0bcb70f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyne\n",
    "import pyne.spatialsolver\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3b9711bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict = {'Name': 'Jane', 'Age': 27}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ef0b0db7",
   "metadata": {},
   "source": [
    "The `spatialsolver` module takes in a dictionary that contains all of the input information required to run the solvers. There are many entries to allow a high degree of customization, not all of which are required.  To find which entries are required, see the spatial solver documentation in the python api."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "b3b4dce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['solver'] = \"AHOTN\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "035d94d0",
   "metadata": {},
   "source": [
    "The neutron transport equations can be solved using various methods, and the spatial solver supports three different families of solvers. These families, along with their corresponding keys, are described in the theory manual as follows:\n",
    "\n",
    "- **AHOTN** - Arbitrarily Higher Order Transport Method: This method involves solving the neutron transport equations using a higher-order scheme, allowing for increased accuracy in the solution.\n",
    "\n",
    "- **DGFEM** - Discontinuous Galerkin Finite Element Method: The DGFEM method employs a finite element approach with discontinuous basis functions to solve the transport equations. This method is particularly useful for handling complex geometries and capturing discontinuities in the solution.\n",
    "\n",
    "- **SCTSTEP** - SCT Step Algorithm: This algorithm, similar to Duo's SCT algorithm, is implemented in three-dimensional Cartesian geometry. It provides a specific approach for solving the neutron transport equations in this geometric configuration.\n",
    "\n",
    "By selecting the appropriate \"solver\" key, you can choose which family of solvers you would like to use for solving the neutron transport equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "f2182a94",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['solver_type'] = \"LN\" "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e9531bc5",
   "metadata": {},
   "source": [
    "Among the families of solvers available, except for SCTSTEP, each family offers multiple choices for solving the neutron transport equation. Here are the specific options for each family:\n",
    "\n",
    "AHOTN:\n",
    "\n",
    "- LN - Arbitrarily Higher Order Transport Method of the Nodal Type, using the linear-nodal method.\n",
    "- LL - Arbitrarily Higher Order Transport Method of the Nodal Type, using the linear-linear method.\n",
    "- NEFD - Arbitrarily Higher Order Transport Method of the Nodal Type, utilizing the unknown nodal flux moments (NEFD algorithm).\n",
    "\n",
    "DGFEM:\n",
    "\n",
    "- LD - Discontinuous Galerkin Finite Element Method (DGFEM) with a linear discontinuous (LD) approximation for the angular flux.\n",
    "- DENSE - Discontinuous Galerkin Finite Element Method (DGFEM) using dense Lagrange polynomials.\n",
    "- LAGRANGE - Discontinuous Galerkin Finite Element Method (DGFEM) employing Lagrange polynomials.\n",
    "\n",
    "SCTSTEP:\n",
    "- SCT Step algorithm, similar to Duo's SCT algorithm, implemented in three-dimensional Cartesian geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "03e51740",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['spatial_order'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19a85b18",
   "metadata": {},
   "source": [
    "The Spatial expansion order is the expansion order of the spatial moment.  It is also known as lambda, and for all AHOTN solvers it must be 0, 1 or 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "722bca70",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['angular_quadrature_order'] = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fe26b06",
   "metadata": {},
   "source": [
    "  The angular quadrature order is the number of angles to be used per octant.  \n",
    "  For N sets of angles, there will be (N * (N + 2) / 8) ordinates per octant. \n",
    "  The quadrature order may only be an even number!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "a3c6a2c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['angular_quadrature_type'] = 1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e89e6be0",
   "metadata": {},
   "source": [
    "  The quadrature type is the type of quadrature scheme the code uses.  \n",
    "  The possibilities are:\n",
    "  \n",
    "- TWOTRAN\n",
    "- EQN\n",
    "- Read-in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "25af56ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['nodes_xyz'] = [4,4,4]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c6420684",
   "metadata": {},
   "source": [
    "The variable `nodes_xyz` represents the number of nodes in the x, y, and z directions. It should be stored as a 1D array with three entries in the following order:\n",
    "\n",
    "- The first entry corresponds to the number of nodes in the x direction (integer).\n",
    "- The second entry corresponds to the number of nodes in the y direction (integer).\n",
    "- The third entry corresponds to the number of nodes in the z direction (integer)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "58931647",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['num_groups'] = 1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0b5c16b4",
   "metadata": {},
   "source": [
    "`num_groups` specifies the number of material groups you are using in the material id and cross section files found in later entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "7aba2e8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['num_materials'] = 1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ab119756",
   "metadata": {},
   "source": [
    "`num_materials` is the number of different materials used in the mesh ('material_id')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "75b80c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['x_cells_widths'] = [0.25, 0.25, 0.25, 0.25]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "5628c01c",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['y_cells_widths'] = [0.25, 0.25, 0.25, 0.25]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "b9bed94b",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['z_cells_widths'] = [0.25, 0.25, 0.25, 0.25]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "45293347",
   "metadata": {},
   "source": [
    "The variables `x_cells_widths`, `y_cells_widths`, and `z_cells_widths` represent the widths of each cell in the x, y, and z directions, respectively. In order to ensure that adjacent cells fit together properly, each unique cell cannot have a unique size. Instead, the cell width specified applies to all cells in the plane perpendicular to the corresponding axis.\n",
    "\n",
    "To clarify, if you assign a value of 1 to the first entry in the `x_cells_widths` array, it means that all cells with an x dimension of 1 will have a width of 1 unit.\n",
    "\n",
    "Please note that the input arrays for `x_cells_widths`, `y_cells_widths`, and `z_cells_widths` should have dimensions of 1 by the number of nodes in the corresponding axis, and all entries in the array should be filled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "32fa6e06",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['x_boundry_conditions'] = [2, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "8b8f79bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['y_boundry_conditions'] = [2, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "9fd8124d",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['z_boundry_conditions'] = [2, 2]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e3b330ae",
   "metadata": {},
   "source": [
    "The variables `x_boundary_conditions`, `y_boundary_conditions`, and `z_boundary_conditions` represent the boundary conditions for each face of a cubic mesh. The specific entries are defined as follows:\n",
    "\n",
    "`x_boundary_conditions`:\n",
    "\n",
    "- 0: Boundary condition for the start face in the x-direction (xsbc).\n",
    "\n",
    "- 1: Boundary condition for the end face in the x-direction (xebc).\n",
    "\n",
    "`y_boundary_conditions`:\n",
    "\n",
    "- 0: Boundary condition for the start face in the y-direction (ysbc).\n",
    "\n",
    "- 1: Boundary condition for the end face in the y-direction (yebc).\n",
    "\n",
    "`z_boundary_conditions`:\n",
    "\n",
    "- 0: Boundary condition for the start face in the z-direction (zsbc).\n",
    "\n",
    "- 1: Boundary condition for the end face in the z-direction (zebc).\n",
    "\n",
    "The supported boundary conditions are as follows:\n",
    "\n",
    "0 - Vacuum: Represents a vacuum boundary condition.\n",
    "\n",
    "1 - Reflective: Indicates a reflective boundary condition.\n",
    "\n",
    "2 - Fixed Inflow: Denotes a fixed inflow boundary condition.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "0547e4e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['material_id'] = [ [ [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1] ], \n",
    "                              [ [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1] ],  \n",
    "                              [ [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1] ],  \n",
    "                              [ [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1] ] ]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1234eaef",
   "metadata": {},
   "source": [
    "The array `material_id` holds information about the materials used in a cubic mesh, which is intended for solving the neutron transport method. Please ensure that the dimensions of the array match the dimensions of the mesh cells, so that each spatial cell corresponds to a single material number. The cells should be ordered in the x, y, and z directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "31db4444",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['quadrature_file'] = 'quad_file'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1624ecef",
   "metadata": {},
   "source": [
    "The `quad_file` refers to the quadrature file used in the neutron transport method. It is necessary only when the quadrature type is set to 2. In that case, providing a valid quadrature file is mandatory. If your quadrature type is different from 2, you can create an empty file and pass it as input for this parameter. Please refer to the formatting notes in the Spatial Solver Python API for further guidance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "4cf40fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['xs_file'] = 'xs'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cec224b1",
   "metadata": {},
   "source": [
    "The `xs_file` refers to the file that contains the cross-sectional data for the materials in your mesh, represented by the `material_id` array. The cross-sectional data should be formatted similar to the following example for a two-material xs file:\n",
    "      \n",
    "    ! Cross section file\n",
    "    ! Material # 1\n",
    "    ! Group #1\n",
    "    5.894     ! Total XS\n",
    "    1.8       ! Scattering matrix\n",
    "    ! Material # 2\n",
    "    ! Group #1\n",
    "    1.237      ! Total XS\n",
    "    0.12       ! Scattering matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "b3e70f9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['source_input_file'] = 'src_4.dat'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e7a7db8b",
   "metadata": {},
   "source": [
    "Note: See input file formatting notes in the Source File Formatting section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "754b7641",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['bc_input_file'] = 'bc_4.dat'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0d4d6719",
   "metadata": {},
   "source": [
    "The `bc_input_file` is the file used to specify the boundary conditions for the neutron inflow on the faces of the mesh. Specifically, it is used for the boundary conditions that are designated as 2 (fixed inflow). To correctly format the boundary condition input file, please refer to the Boundry Condition formatting notes in the Spatial Solver Python API for detailed instructions and guidelines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "081f757d",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['flux_output_file'] = 'phi_4.ahot'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8e11ee0a",
   "metadata": {},
   "source": [
    "`flux_output_file` is the output file for the angular flux to be printed to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "62edf840",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['convergence_criterion'] = 1.e-12"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b029eb21",
   "metadata": {},
   "source": [
    "The calculation is considered complete and the solution is considered converged when the flux in each cell at the current iteration is within a certain threshold called the `convergence_criterion` of the previous iteration. Typically, the convergence criterion is defined as the relative difference between the current and previous iterates. However, in situations where the flux values are very small, the absolute difference is used instead. For more details on the convergence tolerance and its determination, please refer to the Convergence Tolerance entry in the documentation or relevant resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "c39c9a29",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['max_iterations'] = 6000"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "29ea6c7f",
   "metadata": {},
   "source": [
    "`max_iterations` is the maximum number of times the mesh should be sweeped.\n",
    "\n",
    "**Note:** If this number of iterations is reached before the convergence criterion is satisfied, the calculation will terminate and report the current flux estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "bc9da608",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['moments_converged'] = 0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ed5ee0ec",
   "metadata": {},
   "source": [
    "`moments_converged` is the number of moments that should be converged upon for each quadrature in the solution space. Value for moments converged must be in range [0, `spatial_order_in`]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "cb8a63ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict['converge_tolerence'] = 1.e-10"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7d22d0c6",
   "metadata": {},
   "source": [
    "Converge tolerance is the tolerance that determines how the difference between flux iterates (df) that is used to determine convergence will be calculated. \n",
    "\n",
    "df is calculated as follows:\n",
    "\n",
    "```python\n",
    "   f = current flux value\n",
    "   ct = convergence tolerance (value for this key, \"converge_tolerance\")\n",
    "   f1 = flux value from the previous iteration\n",
    "   If f1 > ct:\n",
    "     df = absolute(f - f1) / f1\n",
    "   Else\n",
    "     df = absolute(f - f1)\n",
    "```\n",
    "The idea is to use the absolute difference instead of the relative difference between iterates when the flux is very small to help avoid rounding error.</pre>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "326bd54e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_results = {}\n",
    "dict_results = pyne.spatialsolver.solve(input_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6178e89",
   "metadata": {},
   "source": [
    "Before doing anything with the resulting data, you should check if the solver succesfully ran.\n",
    "If the dictionary key 'success' is 1 (true), the job ran succesfully.  If it is 0 (false), you\n",
    "were not so succesfull."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "a988d02a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Yay, job ran succesfully!\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "if(dict_results['success']):\n",
    "    print('Yay, job ran succesfully!')\n",
    "print(dict_results['success'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b6efb8b",
   "metadata": {},
   "source": [
    "If you were not so lucky, and your job failed, the following key will give you the error message.  It will be 0 if the codes ran succesfully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "dbf76827",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "print(dict_results['error_msg'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea68461d",
   "metadata": {},
   "source": [
    "To get the results of the solver, create a output dictionary to store all the data from the solver, and then use solve to populate it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "766e4663",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[3.5265019965239373, 3.0926025752216213, 3.0926025752216213, 3.5265019965239373], [3.0926025752216204, 2.7320973265457127, 2.7320973265457127, 3.092602575221621], [3.092602575221621, 2.732097326545713, 2.732097326545713, 3.092602575221621], [3.526501996523937, 3.0926025752216213, 3.0926025752216213, 3.5265019965239364]], [[2.890218323639852, 2.612848109027689, 2.6128481090276896, 2.890218323639852], [2.6128481090276896, 2.385716782664984, 2.385716782664984, 2.6128481090276896], [2.6128481090276896, 2.385716782664984, 2.385716782664984, 2.61284810902769], [2.890218323639852, 2.61284810902769, 2.6128481090276896, 2.890218323639852]], [[2.8902183236398526, 2.61284810902769, 2.61284810902769, 2.8902183236398526], [2.6128481090276896, 2.3857167826649843, 2.3857167826649843, 2.6128481090276896], [2.6128481090276896, 2.3857167826649843, 2.3857167826649843, 2.6128481090276896], [2.8902183236398526, 2.61284810902769, 2.61284810902769, 2.8902183236398526]], [[3.526501996523937, 3.0926025752216213, 3.0926025752216213, 3.5265019965239364], [3.092602575221621, 2.732097326545713, 2.7320973265457127, 3.092602575221621], [3.092602575221621, 2.7320973265457127, 2.7320973265457127, 3.092602575221621], [3.526501996523937, 3.0926025752216213, 3.092602575221621, 3.5265019965239373]]]\n"
     ]
    }
   ],
   "source": [
    "print(dict_results['flux'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f109e3c2",
   "metadata": {},
   "source": [
    "To obtain the total time taken by your job, you can retrieve it using the `total_time` key from the dictionary. This key provides the information about the elapsed time for the entire job execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "b87a373a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total solving time was: \n",
      "0.004986000000000157\n"
     ]
    }
   ],
   "source": [
    "print('Total solving time was: ')\n",
    "print(dict_results['total_time'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7ddd88fb",
   "metadata": {},
   "source": [
    "To obtain the absolute system time when the solver call finished, you can calculate it by adding the total job time (retrieved using the `total_time` key) and the solver start time (retrieved using the `time_start` key). The `time_start` key provides the system time when the solver call began, and by adding the total job time to the start time, you can determine the absolute system time when the solver call finished."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "7eb9f2e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver call started at: \n",
      "2.292049\n",
      "Solver call finished at: \n",
      "2.297035\n"
     ]
    }
   ],
   "source": [
    "print('Solver call started at: ')\n",
    "print(dict_results['time_start'])\n",
    "print('Solver call finished at: ')\n",
    "print(dict_results['time_start'] + dict_results['total_time'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
